home *** CD-ROM | disk | FTP | other *** search
/ Aminet 6 / Aminet 6 - June 1995.iso / Aminet / gfx / 3d / irit50src.lha / irit5 / cagd_lib / sbsp_aux.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-16  |  45.3 KB  |  1,271 lines
  1. /******************************************************************************
  2. * SBsp-Aux.c - Bspline surface auxilary routines.                  *
  3. *******************************************************************************
  4. * Written by Gershon Elber, July. 90.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "cagd_loc.h"
  11.  
  12. /* Define some marcos to make some of the routines below look better. They  */
  13. /* calculate the index of the U, V point of the control mesh in Points.        */
  14. #define DERIVED_SRF(U, V)    CAGD_MESH_UV(DerivedSrf, U, V)
  15. #define RAISED_SRF(U, V)    CAGD_MESH_UV(RaisedSrf, U, V)
  16. #define SRF(U, V)        CAGD_MESH_UV(Srf, U, V)
  17.  
  18. #define NORMAL_EPSILON        1e-4
  19.  
  20. static CagdBType EvalTangentVector(CagdVType Vec,
  21.                    CagdCrvStruct *Crv,
  22.                    CagdRType t,
  23.                    CagdRType TMin,
  24.                    CagdRType TMax);
  25.  
  26. /*****************************************************************************
  27. * DESCRIPTION:                                                               M
  28. * Given a Bspline surface - subdivides it into two sub-surfaces at the given M
  29. * parametric value.                                                          M
  30. *   Returns pointer to first surface in a list of two subdivided surfaces.   M
  31. *                                                                            *
  32. * PARAMETERS:                                                                M
  33. *   Srf:      To subdivide at parameter value t.                             M
  34. *   t:        Parameter value to subdivide Srf at.                           M
  35. *   Dir:      Direction of subdivision. Either U or V.                       M
  36. *                                                                            *
  37. * RETURN VALUE:                                                              M
  38. *   CagdSrfStruct *:  A list of the two subdivided surfaces.                 M
  39. *                                                                            *
  40. * KEYWORDS:                                                                  M
  41. *   BspSrfSubdivAtParam, subdivision, refinement                             M
  42. *****************************************************************************/
  43. CagdSrfStruct *BspSrfSubdivAtParam(CagdSrfStruct *Srf,
  44.                    CagdRType t,
  45.                    CagdSrfDirType Dir)
  46. {
  47.     CagdBType
  48.     NewSrf = FALSE,
  49.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Srf);
  50.     int i, j, Row, Col, KVLen, Index1, Index2, Mult,
  51.     LULength, RULength, LVLength, RVLength,    ULength, VLength,
  52.     UOrder = Srf -> UOrder,
  53.     VOrder = Srf -> VOrder,
  54.     MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
  55.     CagdRType *RefKV, **Pts, **LPts, **RPts, *LOnePts, *ROnePts, *OnePts,
  56.         UMin, UMax, VMin, VMax;
  57.     CagdSrfStruct *RSrf, *LSrf;
  58.     BspKnotAlphaCoeffType *A;
  59.  
  60.     if (CAGD_IS_PERIODIC_SRF(Srf)) {
  61.     NewSrf = TRUE;
  62.     Srf = CnvrtPeriodic2FloatSrf(Srf);
  63.     }
  64.  
  65.     ULength = Srf -> ULength;
  66.     VLength = Srf -> VLength;
  67.  
  68.     switch (Dir) {
  69.     case CAGD_CONST_U_DIR:
  70.         RefKV = Srf -> UKnotVector;
  71.         KVLen = UOrder + ULength;
  72.         Index1 = BspKnotLastIndexL(RefKV, KVLen, t);
  73.         if (Index1 + 1 < UOrder)
  74.         Index1 = UOrder - 1;
  75.         Index2 = BspKnotFirstIndexG(RefKV, KVLen, t);
  76.         if (Index2 > ULength)
  77.         Index2 = ULength;
  78.         LSrf = BspSrfNew(Index1 + 1, VLength,
  79.                  UOrder, VOrder, Srf -> PType);
  80.         RSrf = BspSrfNew(ULength - Index2 + UOrder, VLength,
  81.                  UOrder, VOrder, Srf -> PType);
  82.         Mult = UOrder - 1 - (Index2 - Index1 - 1);
  83.  
  84.         /* Update the new knot vectors. */
  85.         CAGD_GEN_COPY(LSrf -> UKnotVector,
  86.               Srf -> UKnotVector,
  87.               sizeof(CagdRType) * (Index1 + 1));
  88.         /* Close the knot vector with multiplicity Order: */
  89.         for (j = Index1 + 1; j <= Index1 + UOrder; j++)
  90.         LSrf -> UKnotVector[j] = t;
  91.         CAGD_GEN_COPY(&RSrf -> UKnotVector[UOrder],
  92.               &Srf -> UKnotVector[Index2],
  93.               sizeof(CagdRType) * (ULength + UOrder - Index2));
  94.         /* Make sure knot vector starts with multiplicity Order: */
  95.         for (j = 0; j < UOrder; j++)
  96.         RSrf -> UKnotVector[j] = t;
  97.         
  98.         /* And copy the other direction knot vectors. */
  99.         CAGD_GEN_COPY(LSrf -> VKnotVector,
  100.               Srf -> VKnotVector,
  101.               sizeof(CagdRType) * (VOrder + VLength));
  102.         CAGD_GEN_COPY(RSrf -> VKnotVector,
  103.               Srf -> VKnotVector,
  104.               sizeof(CagdRType) * (VOrder + VLength));
  105.         break;
  106.     case CAGD_CONST_V_DIR:
  107.         RefKV = Srf -> VKnotVector;
  108.         KVLen = VOrder + VLength;
  109.         Index1 = BspKnotLastIndexL(RefKV, KVLen, t);
  110.         if (Index1 + 1 < VOrder)
  111.         Index1 = VOrder - 1;
  112.         Index2 = BspKnotFirstIndexG(RefKV, KVLen, t);
  113.         if (Index2 > VLength)
  114.         Index2 = VLength;
  115.         LSrf = BspSrfNew(ULength, Index1 + 1,
  116.                  UOrder, VOrder, Srf -> PType);
  117.         RSrf = BspSrfNew(ULength, VLength - Index2 + VOrder,
  118.                  UOrder, VOrder, Srf -> PType);
  119.         Mult = VOrder - 1 - (Index2 - Index1 - 1);
  120.  
  121.         /* Update the new knot vectors. */
  122.         CAGD_GEN_COPY(LSrf -> VKnotVector,
  123.               Srf -> VKnotVector,
  124.               sizeof(CagdRType) * (Index1 + 1));
  125.         /* Close the knot vector with multiplicity Order: */
  126.         for (j = Index1 + 1; j <= Index1 + VOrder; j++)
  127.         LSrf -> VKnotVector[j] = t;
  128.         CAGD_GEN_COPY(&RSrf -> VKnotVector[VOrder],
  129.               &Srf -> VKnotVector[Index2],
  130.               sizeof(CagdRType) * (VLength + VOrder - Index2));
  131.         /* Make sure knot vector starts with multiplicity Order: */
  132.         for (j = 0; j < VOrder; j++)
  133.         RSrf -> VKnotVector[j] = t;
  134.         
  135.         /* And copy the other direction knot vectors. */
  136.         CAGD_GEN_COPY(LSrf -> UKnotVector,
  137.               Srf -> UKnotVector,
  138.               sizeof(CagdRType) * (UOrder + ULength));
  139.         CAGD_GEN_COPY(RSrf -> UKnotVector,
  140.               Srf -> UKnotVector,
  141.               sizeof(CagdRType) * (UOrder + ULength));
  142.         break;
  143.     default:
  144.         Mult = 1;
  145.         RefKV = NULL;
  146.         LSrf = RSrf = NULL;
  147.         CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  148.         break;
  149.     }
  150.  
  151.     Pts = Srf -> Points;
  152.     LPts = LSrf -> Points;
  153.     RPts = RSrf -> Points;
  154.     LULength = LSrf -> ULength,
  155.     RULength = RSrf -> ULength;
  156.     LVLength = LSrf -> VLength,
  157.     RVLength = RSrf -> VLength;
  158.  
  159.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  160.  
  161.     switch (Dir) {
  162.     case CAGD_CONST_U_DIR:
  163.         /* Compute the Alpha refinement matrix. */
  164.         if (Mult > 0) {
  165.         CagdRType
  166.             *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);
  167.  
  168.         CAGD_DOMAIN_T_VERIFY(t, UMin, UMax);
  169.         if (t == UMax)
  170.             t -= CAGD_DOMAIN_EPSILON;
  171.         for (i = 0; i < Mult; i++)
  172.             NewKV[i] = t;
  173.         A = BspKnotEvalAlphaCoefMerge(UOrder, RefKV, ULength,
  174.                                 NewKV, Mult);
  175.         IritFree((VoidPtr) NewKV);
  176.         }
  177.         else
  178.         A = BspKnotEvalAlphaCoefMerge(UOrder, RefKV, ULength, NULL, 0);
  179.  
  180.         /* Now work on the two new surfaces meshes. */
  181.  
  182.         /* Note that Mult can be negative in cases where original       */
  183.         /* multiplicity was order or more and we need to compensate     */
  184.         /* here, since Alpha matrix will be just a unit matrix then.    */
  185.         Mult = Mult >= 0 ? 0 : -Mult;
  186.  
  187.         for (Row = 0; Row < VLength; Row++) {
  188.             /* Blend Srf into LSrf. */
  189.             for (j = IsNotRational; j <= MaxCoord; j++) {
  190.             LOnePts = &LPts[j][Row * LULength];
  191.             OnePts = &Pts[j][Row * ULength];
  192.             for (i = 0; i < LULength; i++, LOnePts++)
  193.                 CAGD_ALPHA_BLEND(A, i, OnePts, Srf -> ULength, LOnePts);
  194.         }
  195.  
  196.             /* Blend Srf into LSrf. */
  197.         for (j = IsNotRational; j <= MaxCoord; j++) {
  198.             ROnePts = &RPts[j][Row * RULength];
  199.             OnePts = &Pts[j][Row * ULength];
  200.             for (i = LSrf -> ULength - 1 + Mult;
  201.              i < LSrf -> ULength + RSrf -> ULength - 1 + Mult;
  202.              i++, ROnePts++)
  203.             CAGD_ALPHA_BLEND(A, i, OnePts, Srf -> ULength, ROnePts);
  204.         }
  205.         }
  206.  
  207.         BspKnotFreeAlphaCoef(A);
  208.         break;
  209.     case CAGD_CONST_V_DIR:
  210.         /* Compute the Alpha refinement matrix. */
  211.         if (Mult > 0) {
  212.         CagdRType
  213.             *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);
  214.  
  215.         CAGD_DOMAIN_T_VERIFY(t, VMin, VMax);
  216.         if (t == VMax)
  217.             t -= CAGD_DOMAIN_EPSILON;
  218.         for (i = 0; i < Mult; i++)
  219.             NewKV[i] = t;
  220.         A = BspKnotEvalAlphaCoefMerge(VOrder, RefKV, VLength,
  221.                           NewKV, Mult);
  222.         IritFree((VoidPtr) NewKV);
  223.         }
  224.         else
  225.         A = BspKnotEvalAlphaCoefMerge(VOrder, RefKV, VLength, NULL, 0);
  226.  
  227.         /* Now work on the two new surfaces meshes. */
  228.  
  229.         /* Note that Mult can be negative in cases where original       */
  230.         /* multiplicity was order or more and we need to compensate     */
  231.         /* here, since Alpha matrix will be just a unit matrix then.    */
  232.         Mult = Mult >= 0 ? 0 : -Mult;
  233.  
  234.         for (Col = 0; Col < ULength; Col++) {
  235.         LULength = LSrf -> ULength;
  236.         RULength = RSrf -> ULength;
  237.  
  238.         /* Blend Srf into LSrf. */
  239.         for (j = IsNotRational; j <= MaxCoord; j++) {
  240.             LOnePts = &LPts[j][Col];
  241.             OnePts = &Pts[j][Col];
  242.             for (i = 0;
  243.                  i < LVLength;
  244.                  i++, LOnePts += LULength)
  245.                 CAGD_ALPHA_BLEND_STEP(A, i, OnePts, Srf -> VLength,
  246.                           LOnePts, ULength);
  247.         }
  248.  
  249.         /* Blend Srf into RSrf. */
  250.         for (j = IsNotRational; j <= MaxCoord; j++) {
  251.             ROnePts = &RPts[j][Col];
  252.             OnePts = &Pts[j][Col];
  253.             for (i = LVLength - 1 + Mult;
  254.              i < LVLength + RVLength - 1 + Mult;
  255.              i++, ROnePts += RULength)
  256.             CAGD_ALPHA_BLEND_STEP(A, i, OnePts, Srf -> VLength,
  257.                           ROnePts, ULength);
  258.         }
  259.         }
  260.  
  261.         BspKnotFreeAlphaCoef(A);
  262.         break;
  263.     default:
  264.         CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  265.         break;
  266.     }
  267.  
  268.     BspKnotMakeRobustKV(RSrf -> UKnotVector,
  269.             RSrf -> UOrder + RSrf -> ULength);
  270.     BspKnotMakeRobustKV(RSrf -> VKnotVector,
  271.             RSrf -> VOrder + RSrf -> VLength);
  272.  
  273.     BspKnotMakeRobustKV(LSrf -> UKnotVector,
  274.             LSrf -> UOrder + LSrf -> ULength);
  275.     BspKnotMakeRobustKV(LSrf -> VKnotVector,
  276.             LSrf -> VOrder + LSrf -> VLength);
  277.  
  278.     LSrf -> Pnext = RSrf;
  279.  
  280.     if (NewSrf)
  281.     CagdSrfFree(Srf);
  282.  
  283.     return LSrf;
  284. }
  285.  
  286. /*****************************************************************************
  287. * DESCRIPTION:                                                               M
  288. *  Inserts n knot, all with the value t in direction Dir. In no case will    M
  289. * the multiplicity of a knot be greater or equal to the curve order.         M
  290. *                                                                            *
  291. * PARAMETERS:                                                                M
  292. *   Srf:        To refine by insertion (upto) n knot of value t.             M
  293. *   Dir:        Direction of refinement. Either U or V.                      M
  294. *   t:          Parameter value of new knot to insert.                       M
  295. *   n:          Maximum number of times t should be inserted.                M
  296. *                                                                            *
  297. * RETURN VALUE:                                                              M
  298. *   CagdSrfStruct *:  Refined Srf with n knots of value t in direction Dir.  M
  299. *                                                                            *
  300. * KEYWORDS:                                                                  M
  301. *   BspSrfKnotInsertNSame, refinement, subdivision                           M
  302. *****************************************************************************/
  303. CagdSrfStruct *BspSrfKnotInsertNSame(CagdSrfStruct *Srf,
  304.                      CagdSrfDirType Dir,
  305.                      CagdRType t,
  306.                      int n)
  307. {
  308.     CagdBType
  309.     NewSrf = FALSE;
  310.     int i, CrntMult, Mult;
  311.     CagdSrfStruct *RefinedSrf;
  312.  
  313.     if (CAGD_IS_PERIODIC_SRF(Srf)) {
  314.         NewSrf = TRUE;
  315.         Srf = CnvrtPeriodic2FloatSrf(Srf);
  316.     }
  317.  
  318.     switch (Dir) {
  319.     case CAGD_CONST_U_DIR:
  320.         CrntMult = BspKnotFindMult(Srf -> UKnotVector, Srf -> UOrder,
  321.                              Srf -> ULength, t),
  322.         Mult = MIN(n, Srf -> UOrder - CrntMult - 1);
  323.         break;
  324.     case CAGD_CONST_V_DIR:
  325.         CrntMult = BspKnotFindMult(Srf -> VKnotVector, Srf -> VOrder,
  326.                              Srf -> VLength, t),
  327.         Mult = MIN(n, Srf -> VOrder - CrntMult - 1);
  328.         break;
  329.     default:
  330.         Mult = 0;
  331.         CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  332.         break;
  333.     }
  334.  
  335.     if (Mult > 0) {
  336.     CagdRType
  337.         *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);
  338.  
  339.     for (i = 0; i < Mult; i++)
  340.         NewKV[i] = t;
  341.  
  342.     RefinedSrf = BspSrfKnotInsertNDiff(Srf, Dir, FALSE, NewKV, Mult);
  343.  
  344.     IritFree((VoidPtr) NewKV);
  345.     }
  346.     else {
  347.     RefinedSrf = CagdSrfCopy(Srf);
  348.     }
  349.  
  350.     if (NewSrf)
  351.     CagdSrfFree(Srf);
  352.  
  353.     return RefinedSrf;
  354. }
  355.  
  356. /*****************************************************************************
  357. * DESCRIPTION:                                                               M
  358. *  Inserts n knot with different values as defined by the vector t. If,      M
  359. * however, Replace is TRUE, the knot are simply replacing the current knot   M
  360. * vector.                                                                    M
  361. *                                                                            *
  362. * PARAMETERS:                                                                M
  363. *   Srf:        To refine by insertion (upto) n knot of value t.             M
  364. *   Dir:        Direction of refinement. Either U or V.                      M
  365. *   Replace:    if TRUE, the n knots in t should replace the knot vector     M
  366. *               of size n of Srf. Sizes must match. If False, n new knots    M
  367. *               as defined by t will be introduced into Srf.                 M
  368. *   t:          New knots to introduce/replace knot vector of Srf.           M
  369. *   n:          Size of t.                                                   M
  370. *                                                                            *
  371. * RETURN VALUE:                                                              M
  372. *   CagdSrfStruct *:    Refined Srf with n new knots in direction Dir.       M
  373. *                                                                            *
  374. * KEYWORDS:                                                                  M
  375. *   BspSrfKnotInsertNDiff, refinement, subdivision                           M
  376. *****************************************************************************/
  377. CagdSrfStruct *BspSrfKnotInsertNDiff(CagdSrfStruct *Srf,
  378.                      CagdSrfDirType Dir,
  379.                      int Replace,
  380.                      CagdRType *t,
  381.                      int n)
  382. {
  383.     CagdBType
  384.     NewSrf = FALSE,
  385.     IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
  386.     int i, Row, Col, ULength, VLength,
  387.     UOrder = Srf -> UOrder,
  388.     VOrder = Srf -> VOrder,
  389.     MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
  390.     CagdSrfStruct
  391.     *RefSrf = NULL;
  392.  
  393.     if (CAGD_IS_PERIODIC_SRF(Srf)) {
  394.         NewSrf = TRUE;
  395.         Srf = CnvrtPeriodic2FloatSrf(Srf);
  396.     }
  397.  
  398.     ULength = Srf -> ULength;
  399.     VLength = Srf -> VLength;
  400.  
  401.     if (Replace) {
  402.     for (i = 1; i < n; i++)
  403.         if (t[i] < t[i - 1])
  404.         CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);
  405.  
  406.         switch (Dir) {
  407.         case CAGD_CONST_U_DIR:
  408.         if (Srf -> UOrder + Srf -> ULength != n)
  409.             CAGD_FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);
  410.  
  411.         RefSrf = CagdSrfCopy(Srf);
  412.         for (i = 0; i < n; i++)
  413.             RefSrf -> UKnotVector[i] = *t++;
  414.         break;
  415.         case CAGD_CONST_V_DIR:
  416.         if (Srf -> VOrder + Srf -> VLength != n)
  417.             CAGD_FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);
  418.  
  419.         RefSrf = CagdSrfCopy(Srf);
  420.         for (i = 0; i < n; i++)
  421.             RefSrf -> VKnotVector[i] = *t++;
  422.         break;
  423.         default:
  424.         CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  425.         break;
  426.     }
  427.     }
  428.     else if (n == 0) {
  429.     RefSrf = CagdSrfCopy(Srf);
  430.     }
  431.     else {
  432.     int j, LengthKVt, RULength, RVLength;
  433.     BspKnotAlphaCoeffType *A;
  434.     CagdRType *MergedKVt, UMin, UMax, VMin, VMax,
  435.         *UKnotVector = Srf -> UKnotVector,
  436.         *VKnotVector = Srf -> VKnotVector;
  437.  
  438.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  439.  
  440.     for (i = 1; i < n; i++)
  441.         if (t[i] < t[i - 1])
  442.         CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);
  443.  
  444.     switch (Dir) {
  445.         case CAGD_CONST_U_DIR:
  446.             for (i = 0; i < n; i++) {
  447.             CAGD_DOMAIN_T_VERIFY(t[i], UMin, UMax);
  448.             if (t[i] == UMax)
  449.             t[i] -= CAGD_DOMAIN_EPSILON;
  450.         }
  451.  
  452.         /* Compute the Alpha refinement matrix. */
  453.         MergedKVt = BspKnotMergeTwo(UKnotVector, ULength + UOrder,
  454.                         t, n, 0, &LengthKVt);
  455.         A = BspKnotEvalAlphaCoef(UOrder, UKnotVector, ULength,
  456.                      MergedKVt, LengthKVt - UOrder);
  457.  
  458.             RefSrf = BspSrfNew(ULength + n, VLength, UOrder, VOrder,
  459.                                   Srf -> PType);
  460.         IritFree((VoidPtr) RefSrf -> UKnotVector);
  461.         IritFree((VoidPtr) RefSrf -> VKnotVector);
  462.         RefSrf -> UKnotVector = MergedKVt;
  463.         RefSrf -> VKnotVector = BspKnotCopy(Srf -> VKnotVector,
  464.                     Srf -> VLength + Srf -> VOrder);
  465.  
  466.         RULength = RefSrf -> ULength;
  467.  
  468.         /* Update the control mesh */
  469.         for (Row = 0; Row < VLength; Row++) {
  470.             for (j = IsNotRational; j <= MaxCoord; j++) {
  471.             CagdRType
  472.                 *ROnePts = &RefSrf -> Points[j][Row * RULength],
  473.                 *OnePts = &Srf -> Points[j][Row * ULength];
  474.  
  475.             for (i = 0; i < RULength; i++, ROnePts++)
  476.                 CAGD_ALPHA_BLEND(A, i, OnePts, Srf -> ULength,
  477.                          ROnePts);
  478.             }
  479.         }
  480.  
  481.         BspKnotFreeAlphaCoef(A);
  482.         break;
  483.         case CAGD_CONST_V_DIR:
  484.         for (i = 0; i < n; i++) {
  485.             CAGD_DOMAIN_T_VERIFY(t[i], VMin, VMax);
  486.             if (t[i] == VMax)
  487.             t[i] -= CAGD_DOMAIN_EPSILON;
  488.         }
  489.  
  490.         /* Compute the Alpha refinement matrix. */
  491.         MergedKVt = BspKnotMergeTwo(VKnotVector, VLength + VOrder,
  492.                         t, n, 0, &LengthKVt);
  493.         A = BspKnotEvalAlphaCoef(VOrder, VKnotVector, VLength,
  494.                      MergedKVt, LengthKVt - VOrder);
  495.  
  496.             RefSrf = BspSrfNew(ULength, VLength + n, UOrder, VOrder,
  497.                                   Srf -> PType);
  498.         IritFree((VoidPtr) RefSrf -> UKnotVector);
  499.         IritFree((VoidPtr) RefSrf -> VKnotVector);
  500.         RefSrf -> UKnotVector = BspKnotCopy(Srf -> UKnotVector,
  501.                     Srf -> ULength + Srf -> UOrder);
  502.         RefSrf -> VKnotVector = MergedKVt;
  503.  
  504.         RULength = RefSrf -> ULength;
  505.         RVLength = RefSrf -> VLength;
  506.  
  507.         /* Update the control mesh */
  508.         for (Col = 0; Col < ULength; Col++) {
  509.             for (j = IsNotRational; j <= MaxCoord; j++) {
  510.             CagdRType
  511.                 *ROnePts = &RefSrf -> Points[j][Col],
  512.                 *OnePts = &Srf -> Points[j][Col];
  513.  
  514.             for (i = 0; i < RVLength; i++, ROnePts += RULength)
  515.                 CAGD_ALPHA_BLEND_STEP(A, i, OnePts, Srf -> VLength,
  516.                           ROnePts, ULength);
  517.             }
  518.         }
  519.  
  520.         BspKnotFreeAlphaCoef(A);
  521.         break;
  522.         default:
  523.         CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  524.         break;
  525.     }
  526.     }
  527.  
  528.     BspKnotMakeRobustKV(RefSrf -> UKnotVector,
  529.             RefSrf -> UOrder + RefSrf -> ULength);
  530.     BspKnotMakeRobustKV(RefSrf -> VKnotVector,
  531.             RefSrf -> VOrder + RefSrf -> VLength);
  532.  
  533.     if (NewSrf)
  534.     CagdSrfFree(Srf);
  535.  
  536.     return RefSrf;
  537. }
  538.  
  539. /*****************************************************************************
  540. * DESCRIPTION:                                                               M
  541. * Returns a new Bspline surface, identical to the original but with one      M
  542. * degree higher, in the requested direction Dir.                             M
  543. *                                                                            *
  544. * PARAMETERS:                                                                M
  545. *   Srf:        To raise it degree by one.                                   M
  546. *   Dir:        Direction of degree raising. Either U or V.                  M
  547. *                                                                            *
  548. * RETURN VALUE:                                                              M
  549. *   CagdSrfStruct *:  A surface with one degree higher in direction Dir,     M
  550. *                     representing the same geometry as Srf.             M
  551. *                                                                            *
  552. * KEYWORDS:                                                                  M
  553. *   BspSrfDegreeRaise, degree raising                                        M
  554. *****************************************************************************/
  555. CagdSrfStruct *BspSrfDegreeRaise(CagdSrfStruct *Srf, CagdSrfDirType Dir)
  556. {
  557.     CagdBType
  558.     NewSrf = FALSE,
  559.     IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
  560.     int i, i2, j, RaisedLen, Row, Col, Length,
  561.     Order = Dir == CAGD_CONST_V_DIR ? Srf -> UOrder : Srf -> VOrder,
  562.     MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
  563.     CagdSrfStruct
  564.     *RaisedSrf = NULL;
  565.  
  566.     if (CAGD_IS_PERIODIC_SRF(Srf)) {
  567.         NewSrf = TRUE;
  568.         Srf = CnvrtPeriodic2FloatSrf(Srf);
  569.     }
  570.  
  571.     Length = Dir == CAGD_CONST_V_DIR ? Srf -> ULength : Srf -> VLength;
  572.  
  573.     if (Order > 2) {
  574.     CagdSrfStruct *UnitSrf;
  575.     int UKvLen1 = Srf -> UOrder + Srf -> ULength - 1,
  576.         VKvLen1 = Srf -> VOrder + Srf -> VLength - 1;
  577.     CagdRType
  578.         *UKv = Srf -> UKnotVector,
  579.         *VKv = Srf -> VKnotVector;
  580.  
  581.     /* Degree raise by multiplying by a constant 1 linear surface in the */
  582.     /* raised direction and constant 1 constant surface in the other.    */
  583.  
  584.     switch (Dir) {
  585.         case CAGD_CONST_U_DIR:
  586.         UnitSrf = BspSrfNew(1, 2, 1, 2,
  587.                     CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
  588.         for (i = 0; i < 2; i++)
  589.             UnitSrf -> UKnotVector[i] = i > 0 ? UKv[UKvLen1] : UKv[0];
  590.         for (i = 0; i < 4; i++)
  591.             UnitSrf -> VKnotVector[i] = i > 1 ? VKv[VKvLen1] : VKv[0];
  592.         break;
  593.         case CAGD_CONST_V_DIR:
  594.         UnitSrf = BspSrfNew(2, 1, 2, 1,
  595.                     CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
  596.         for (i = 0; i < 4; i++)
  597.             UnitSrf -> UKnotVector[i] = i > 1 ? UKv[UKvLen1] : UKv[0];
  598.         for (i = 0; i < 2; i++)
  599.             UnitSrf -> VKnotVector[i] = i > 0 ? VKv[VKvLen1] : VKv[0];
  600.         break;
  601.         default:
  602.         UnitSrf = NULL;
  603.         CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  604.         break;
  605.     }
  606.     for (i = 1; i <= MaxCoord; i++)
  607.         UnitSrf -> Points[i][0] = UnitSrf -> Points[i][1] = 1.0;
  608.  
  609.     RaisedSrf = BspSrfMult(Srf, UnitSrf);
  610.  
  611.     CagdSrfFree(UnitSrf);
  612.  
  613.     if (NewSrf)
  614.         CagdSrfFree(Srf);
  615.  
  616.     return RaisedSrf;
  617.     }
  618.  
  619.     /* If surface is linear, degree raising means basically to increase the  */
  620.     /* knot multiplicity of each segment by one and add a middle point for   */
  621.     /* each such segment.                             */
  622.     RaisedLen = Length * 2 - 1;
  623.  
  624.     switch (Dir) {
  625.     case CAGD_CONST_U_DIR:
  626.         RaisedSrf = BspSrfNew(Srf -> ULength, RaisedLen,
  627.                   Srf -> UOrder, Order + 1, Srf -> PType);
  628.  
  629.         /* Update the knot vectors. */
  630.         CAGD_GEN_COPY(RaisedSrf -> UKnotVector,
  631.               Srf -> UKnotVector,
  632.               sizeof(CagdRType) * (Srf -> ULength + Srf -> UOrder));
  633.         for (i = 0; i < 3; i++)
  634.         RaisedSrf -> VKnotVector[i] = Srf -> VKnotVector[0];
  635.         for (i = 2, j = 3; i < Length; i++, j += 2)
  636.         RaisedSrf -> VKnotVector[j] = RaisedSrf -> VKnotVector[j + 1] = 
  637.             Srf -> VKnotVector[i];
  638.         for (i = j; i < j + 3; i++)
  639.         RaisedSrf -> VKnotVector[i] = Srf -> VKnotVector[Length];
  640.  
  641.         /* Update the mesh. */
  642.                for (Col = 0; Col < Srf -> ULength; Col++) {
  643.         for (j = IsNotRational; j <= MaxCoord; j++)  /* First point. */
  644.             RaisedSrf -> Points[j][RAISED_SRF(Col, 0)] =
  645.             Srf -> Points[j][SRF(Col, 0)];
  646.  
  647.         for (i = 1, i2 = 1; i < Length; i++, i2 += 2)
  648.             for (j = IsNotRational; j <= MaxCoord; j++) {
  649.             RaisedSrf -> Points[j][RAISED_SRF(Col, i2)] =
  650.                 Srf -> Points[j][SRF(Col, i - 1)] * 0.5 +
  651.                 Srf -> Points[j][SRF(Col, i)] * 0.5;
  652.             RaisedSrf -> Points[j][RAISED_SRF(Col, i2 + 1)] =
  653.                 Srf -> Points[j][SRF(Col, i)];
  654.             }
  655.         }
  656.         break;
  657.     case CAGD_CONST_V_DIR:
  658.         RaisedSrf = BspSrfNew(RaisedLen, Srf -> VLength,
  659.                   Order + 1, Srf -> VOrder, Srf -> PType);
  660.  
  661.         /* Update the knot vectors. */
  662.         CAGD_GEN_COPY(RaisedSrf -> VKnotVector,
  663.               Srf -> VKnotVector,
  664.               sizeof(CagdRType) * (Srf -> VLength + Srf -> VOrder));
  665.         for (i = 0; i < 3; i++)
  666.         RaisedSrf -> UKnotVector[i] = Srf -> UKnotVector[0];
  667.         for (i = 2, j = 3; i < Length; i++, j += 2)
  668.         RaisedSrf -> UKnotVector[j] = RaisedSrf -> UKnotVector[j + 1] = 
  669.             Srf -> UKnotVector[i];
  670.         for (i = j; i < j + 3; i++)
  671.         RaisedSrf -> UKnotVector[i] = Srf -> UKnotVector[Length];
  672.  
  673.         /* Update the mesh. */
  674.                for (Row = 0; Row < Srf -> VLength; Row++) {
  675.         for (j = IsNotRational; j <= MaxCoord; j++)  /* First point. */
  676.             RaisedSrf -> Points[j][RAISED_SRF(0, Row)] =
  677.             Srf -> Points[j][SRF(0, Row)];
  678.  
  679.         for (i = 1, i2 = 1; i < Length; i++, i2 += 2)
  680.             for (j = IsNotRational; j <= MaxCoord; j++) {
  681.             RaisedSrf -> Points[j][RAISED_SRF(i2, Row)] =
  682.                 Srf -> Points[j][SRF(i - 1, Row)] * 0.5 +
  683.                 Srf -> Points[j][SRF(i, Row)] * 0.5;
  684.             RaisedSrf -> Points[j][RAISED_SRF(i2 + 1, Row)] =
  685.                 Srf -> Points[j][SRF(i, Row)];
  686.             }
  687.         }
  688.         break;
  689.     default:
  690.         CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  691.         break;
  692.     }
  693.  
  694.     if (NewSrf)
  695.     CagdSrfFree(Srf);
  696.  
  697.     return RaisedSrf;
  698. }
  699.  
  700. /*****************************************************************************
  701. * DESCRIPTION:                                                               M
  702. * Returns a new surface equal to the given surface, differentiated once in   M
  703. * the direction Dir.                                                         M
  704. *   Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then: M
  705. * Q(i) = (k - 1) * (P(i+1) - P(i)) / (Kv(i + k) - Kv(i + 1)), i = 0 to k-2.  V
  706. * This is applied to all rows/cols of the surface.                 M
  707. *                                                                            *
  708. * PARAMETERS:                                                                M
  709. *   Srf:        To differentiate.                                            M
  710. *   Dir:        Direction of differentiation. Either U or V.                 M
  711. *                                                                            *
  712. * RETURN VALUE:                                                              M
  713. *   CagdSrfStruct *:   Differentiated surface.                               M
  714. *                                                                            *
  715. * KEYWORDS:                                                                  M
  716. *   BspSrfDerive, derivatives, partial derivatives                           M
  717. *****************************************************************************/
  718. CagdSrfStruct *BspSrfDerive(CagdSrfStruct *Srf,
  719.                 CagdSrfDirType Dir)
  720. {
  721.     CagdBType
  722.     NewSrf = FALSE,
  723.     IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
  724.     int i, j, Row, Col, NewUOrder, NewVOrder, NewULength, NewVLength,
  725.     ULength, VLength,
  726.     UOrder = Srf -> UOrder,
  727.     VOrder = Srf -> VOrder,
  728.     MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
  729.     CagdRType **DPoints, *UKv, *VKv, **Points;
  730.     CagdSrfStruct
  731.         *DerivedSrf = NULL;
  732.  
  733.     if (CAGD_IS_PERIODIC_SRF(Srf)) {
  734.         NewSrf = TRUE;
  735.         Srf = CnvrtPeriodic2FloatSrf(Srf);
  736.     }
  737.  
  738.     if (!IsNotRational)
  739.     return BspSrfDeriveRational(Srf, Dir);
  740.  
  741.     ULength = Srf -> ULength;
  742.     VLength = Srf -> VLength;
  743.     UKv = Srf -> UKnotVector;
  744.     VKv = Srf -> VKnotVector;
  745.     Points = Srf -> Points;
  746.  
  747.     switch (Dir) {
  748.     case CAGD_CONST_U_DIR:
  749.         NewULength = UOrder < 2 ? ULength : ULength - 1;
  750.         NewUOrder = MAX(UOrder - 1, 1);
  751.         DerivedSrf = BspSrfNew(NewULength, VLength,
  752.                    NewUOrder, VOrder, Srf -> PType);
  753.         CAGD_GEN_COPY(DerivedSrf -> UKnotVector, &UKv[UOrder < 2 ? 0 : 1],
  754.               sizeof(CagdRType) * (NewULength + NewUOrder));
  755.         CAGD_GEN_COPY(DerivedSrf -> VKnotVector, VKv,
  756.               sizeof(CagdRType) * (VLength + VOrder));
  757.         DPoints = DerivedSrf -> Points;
  758.  
  759.                for (Row = 0; Row < VLength; Row++)
  760.         for (i = 0; i < NewULength; i++) {
  761.             CagdRType
  762.             Denom = UKv[i + UOrder] - UKv[i + 1];
  763.  
  764.             for (j = IsNotRational; j <= MaxCoord; j++) {
  765.             if (APX_EQ(Denom, 0.0))
  766.                 Denom = INFINITY;
  767.  
  768.             DPoints[j][DERIVED_SRF(i, Row)] =
  769.                 UOrder < 2 ? 0.0
  770.                        : (UOrder - 1) *
  771.                         (Points[j][SRF(i + 1, Row)] -
  772.                          Points[j][SRF(i, Row)]) / Denom;
  773.             }
  774.         }
  775.         break;
  776.     case CAGD_CONST_V_DIR:
  777.         NewVLength = VOrder < 2 ? VLength : VLength - 1;
  778.         NewVOrder = MAX(VOrder - 1, 1);
  779.         DerivedSrf = BspSrfNew(ULength, NewVLength,
  780.                        UOrder, NewVOrder, Srf -> PType);
  781.         CAGD_GEN_COPY(DerivedSrf -> UKnotVector, UKv,
  782.               sizeof(CagdRType) * (ULength + UOrder));
  783.         CAGD_GEN_COPY(DerivedSrf -> VKnotVector, &VKv[VOrder < 2 ? 0 : 1],
  784.               sizeof(CagdRType) * (NewVLength + NewVOrder));
  785.         DPoints = DerivedSrf -> Points;
  786.  
  787.         for (Col = 0; Col < ULength; Col++)
  788.         for (i = 0; i < NewVLength; i++) {
  789.             CagdRType
  790.             Denom = VKv[i + VOrder] - VKv[i + 1];
  791.  
  792.             for (j = IsNotRational; j <= MaxCoord; j++) {
  793.             if (APX_EQ(Denom, 0.0))
  794.                 Denom = INFINITY;
  795.  
  796.             DPoints[j][DERIVED_SRF(Col, i)] =
  797.                 VOrder < 2 ? 0.0
  798.                        : (VOrder - 1) *
  799.                         (Points[j][SRF(Col, i + 1)] -
  800.                          Points[j][SRF(Col, i)]) / Denom;
  801.             }
  802.         }
  803.         break;
  804.     default:
  805.         CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  806.         break;
  807.     }
  808.  
  809.     if (NewSrf)
  810.     CagdSrfFree(Srf);
  811.  
  812.     return DerivedSrf;
  813. }
  814.  
  815. /*****************************************************************************
  816. * DESCRIPTION:                                                               M
  817. * Evaluates the unit tangent to a surface at a given parametric location     M
  818. ( u, v) and given direction Dir.                                             M
  819. *                                                                            *
  820. * PARAMETERS:                                                                M
  821. *   Srf:        Bspline surface to evaluate tangent vector for.              M
  822. *   u, v:       Parametric location of required unit tangent.                M
  823. *   Dir:        Direction of tangent vector. Either U or V.                  M
  824. *                                                                            *
  825. * RETURN VALUE:                                                              M
  826. *   CagdVecStruct *:  A pointer to a static vector holding the unit tangent  M
  827. *                     information.                                           M
  828. *                                                                            *
  829. * KEYWORDS:                                                                  M
  830. *   BspSrfTangent, tangent                                                   M
  831. *****************************************************************************/
  832. CagdVecStruct *BspSrfTangent(CagdSrfStruct *Srf,
  833.                  CagdRType u,
  834.                  CagdRType v,
  835.                  CagdSrfDirType Dir)
  836. {
  837.     CagdVecStruct
  838.     *Tangent = NULL;
  839.     CagdCrvStruct *Crv;
  840.  
  841.     switch (Dir) {
  842.     case CAGD_CONST_V_DIR:
  843.         Crv = BspSrfCrvFromSrf(Srf, v, Dir);
  844.         Tangent = BspCrvTangent(Crv, u);
  845.         CagdCrvFree(Crv);
  846.         break;
  847.     case CAGD_CONST_U_DIR:
  848.         Crv = BspSrfCrvFromSrf(Srf, u, Dir);
  849.         Tangent = BspCrvTangent(Crv, v);
  850.         CagdCrvFree(Crv);
  851.         break;
  852.     default:
  853.         CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
  854.         break;
  855.     }
  856.  
  857.     return Tangent;
  858. }
  859.  
  860. /*****************************************************************************
  861. * DESCRIPTION:                                                               M
  862. * Evaluate the unit normal of a surface at a given parametric location.      M
  863. *   If we fail to compute the normal at given location we retry by moving a  M
  864. * tad.                                                                       M
  865. *                                                                            *
  866. * PARAMETERS:                                                                M
  867. *   Srf:        Bspline surface to evaluate normal vector for.               M
  868. *   u, v:       Parametric location of required unit normal.                 M
  869. *                                                                            *
  870. * RETURN VALUE:                                                              M
  871. *   CagdVecStruct *:  A pointer to a static vector holding the unit normal   M
  872. *                     information.                                           M
  873. *                                                                            *
  874. * KEYWORDS:                                                                  M
  875. *   BspSrfNormal, normal                                                     M
  876. *****************************************************************************/
  877. CagdVecStruct *BspSrfNormal(CagdSrfStruct *Srf, CagdRType u, CagdRType v)
  878. {
  879.     static CagdVecStruct Normal;
  880.     CagdVecStruct *V, T1, T2;
  881.     CagdRType UMin, UMax, VMin, VMax;
  882.  
  883.     CAGD_DOMAIN_GET_AND_VERIFY_SRF(u, v, Srf, UMin, UMax, VMin, VMax);
  884.  
  885.     V = BspSrfTangent(Srf, u, v, CAGD_CONST_U_DIR);
  886.     if (CAGD_LEN_VECTOR(*V) < IRIT_EPSILON)
  887.     V = BspSrfTangent(Srf,
  888.               u > UMin + EPSILON ? u - EPSILON : u + EPSILON,
  889.               v > VMin + EPSILON ? v - EPSILON : v + EPSILON,
  890.               CAGD_CONST_U_DIR);
  891.     CAGD_COPY_VECTOR(T1, *V);
  892.  
  893.     V = BspSrfTangent(Srf, u, v, CAGD_CONST_V_DIR);
  894.     if (CAGD_LEN_VECTOR(*V) < IRIT_EPSILON)
  895.     V = BspSrfTangent(Srf,
  896.               u > UMin + EPSILON ? u - EPSILON : u + EPSILON,
  897.               v > VMin + EPSILON ? v - EPSILON : v + EPSILON,
  898.               CAGD_CONST_V_DIR);
  899.     CAGD_COPY_VECTOR(T2, *V);
  900.  
  901.     /* The normal is the cross product of T1 and T2: */
  902.     Normal.Vec[0] = T1.Vec[1] * T2.Vec[2] - T1.Vec[2] * T2.Vec[1];
  903.     Normal.Vec[1] = T1.Vec[2] * T2.Vec[0] - T1.Vec[0] * T2.Vec[2];
  904.     Normal.Vec[2] = T1.Vec[0] * T2.Vec[1] - T1.Vec[1] * T2.Vec[0];
  905.  
  906.     CAGD_NORMALIZE_VECTOR(Normal);            /* Normalize the vector. */
  907.  
  908.     return &Normal;
  909. }
  910.  
  911. /*****************************************************************************
  912. * DESCRIPTION:                                                               M
  913. * Evaluates the unit normals of a surface at a mesh defined by subdividing   M
  914. * the parametric space into a grid of size UFineNess by VFineNess.         M
  915. *   The normals are saved in a linear CagdVecStruct vector which is          M
  916. * allocated dynamically. Data is saved u inc. first.                 M
  917. *   This routine is much faster than evaluating normal for each point,       M
  918. * individually.                                                              M
  919. *                                                                            *
  920. * PARAMETERS:                                                                M
  921. *   Srf:          To compute normals on a grid of its parametric domain.     M
  922. *   UFineNess:    U Fineness of imposed grid on Srf's parametric domain.     M
  923. *   VFineNess:    V Fineness of imposed grid on Srf's parametric domain.     M
  924. *                                                                            *
  925. * RETURN VALUE:                                                              M
  926. *   CagdVecStruct *:   An vector of unit normals (u increments first).       M
  927. *                                                                            *
  928. * KEYWORDS:                                                                  M
  929. *   BspSrfMeshNormals, normal                                                M
  930. *****************************************************************************/
  931. CagdVecStruct *BspSrfMeshNormals(CagdSrfStruct *Srf,
  932.                  int UFineNess,
  933.                  int VFineNess)
  934. {
  935.     int i, j;
  936.     CagdRType UMin, UMax, VMin, VMax;
  937.     CagdVecStruct *Normals, *NPtr;
  938.     CagdSrfStruct
  939.     *DuSrf = CagdSrfDerive(Srf, CAGD_CONST_U_DIR),
  940.     *DvSrf = CagdSrfDerive(Srf, CAGD_CONST_V_DIR);
  941.  
  942.     UFineNess = MAX(2, UFineNess);
  943.     VFineNess = MAX(2, VFineNess);
  944.  
  945.     Normals = CagdVecArrayNew(UFineNess * VFineNess);
  946.  
  947.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  948.  
  949.     NPtr = Normals;
  950.     for (i = 0; i < UFineNess; i++) {
  951.     CagdRType
  952.        U = UMin + (UMax - UMin) * i / (UFineNess - 1);
  953.     CagdCrvStruct
  954.        *DuCrv2 = NULL,
  955.        *DvCrv2 = NULL,
  956.        *DuCrv = CagdCrvFromSrf(DuSrf, U, CAGD_CONST_U_DIR),
  957.        *DvCrv = CagdCrvFromSrf(DvSrf, U, CAGD_CONST_U_DIR);
  958.  
  959.     for (j = 0; j < VFineNess; j++) {
  960.         CagdVType Du, Dv;
  961.         CagdRType
  962.         V = VMin + (VMax - VMin) * j / (VFineNess - 1);
  963.  
  964.         if (!EvalTangentVector(Du, DuCrv, V, VMin, VMax)) {
  965.         if (DuCrv2 == NULL) {
  966.             DuCrv2 = CagdCrvFromSrf(DuSrf,
  967.                         U > UMin + NORMAL_EPSILON ?
  968.                             U - NORMAL_EPSILON :
  969.                             U + NORMAL_EPSILON,
  970.                         CAGD_CONST_U_DIR);
  971.         }
  972.         EvalTangentVector(Du, DuCrv2, V, VMin, VMax);
  973.         }
  974.  
  975.         if (!EvalTangentVector(Dv, DvCrv, V, VMin, VMax)) {
  976.         if (DvCrv2 == NULL) {
  977.             DvCrv2 = CagdCrvFromSrf(DvSrf,
  978.                         U > UMin + NORMAL_EPSILON ?
  979.                             U - NORMAL_EPSILON :
  980.                             U + NORMAL_EPSILON,
  981.                         CAGD_CONST_U_DIR);
  982.         }
  983.         EvalTangentVector(Dv, DvCrv2, V, VMin, VMax);
  984.         }
  985.  
  986.         CROSS_PROD(NPtr -> Vec, Dv, Du);
  987.         NPtr++;
  988.     }
  989.  
  990.     if (DuCrv2 != NULL)
  991.         CagdCrvFree(DuCrv2);
  992.     if (DvCrv2 != NULL)
  993.         CagdCrvFree(DvCrv2);
  994.     CagdCrvFree(DuCrv);
  995.     CagdCrvFree(DvCrv);
  996.     }
  997.  
  998.     /* Normalize the results. */
  999.     NPtr = Normals;
  1000.     for (i = 0; i < UFineNess; i++) {
  1001.     for (j = 0; j < VFineNess; j++) {
  1002.         CagdRType
  1003.         Len = PT_LENGTH(NPtr -> Vec);
  1004.  
  1005.         if (Len < PT_NORMALIZE_ZERO) {
  1006.         CagdVecStruct *NPtrTmp;
  1007.  
  1008.         /* Try its neighbors. */
  1009.         if (i > 0) {
  1010.             NPtrTmp = NPtr - VFineNess;
  1011.             PT_COPY(NPtr -> Vec, NPtrTmp -> Vec);
  1012.             Len = PT_LENGTH(NPtr -> Vec);
  1013.         }
  1014.         if (Len < PT_NORMALIZE_ZERO && i < UFineNess - 1) {
  1015.             NPtrTmp = NPtr + VFineNess;
  1016.             PT_COPY(NPtr -> Vec, NPtrTmp -> Vec);
  1017.             Len = PT_LENGTH(NPtr -> Vec);
  1018.         }
  1019.         if (Len < PT_NORMALIZE_ZERO && j > 0) {
  1020.             NPtrTmp = NPtr - 1;
  1021.             PT_COPY(NPtr -> Vec, NPtrTmp -> Vec);
  1022.             Len = PT_LENGTH(NPtr -> Vec);
  1023.         }
  1024.         if (Len < PT_NORMALIZE_ZERO && j < VFineNess - 1) {
  1025.             NPtrTmp = NPtr + 1;
  1026.             PT_COPY(NPtr -> Vec, NPtrTmp -> Vec);
  1027.             Len = PT_LENGTH(NPtr -> Vec);
  1028.         }
  1029.         if (Len > PT_NORMALIZE_ZERO) {
  1030.             Len = 1 / Len;
  1031.             PT_SCALE(NPtr -> Vec, Len);
  1032.         }
  1033.         else {
  1034.             /* Do something. */
  1035.             NPtr -> Vec[0] = NPtr -> Vec[1] = 0.0;
  1036.             NPtr -> Vec[2] = 1.0;
  1037.         }
  1038.         }
  1039.         else {
  1040.         Len = 1 / Len;
  1041.         PT_SCALE(NPtr -> Vec, Len);
  1042.         }
  1043.         NPtr++;
  1044.     }
  1045.     }
  1046.  
  1047.  
  1048.     CagdSrfFree(DuSrf);
  1049.     CagdSrfFree(DvSrf);
  1050.  
  1051.     return Normals;
  1052. }
  1053.  
  1054. /*****************************************************************************
  1055. * DESCRIPTION:                                                               *
  1056. *   Evaluates a vector field Crv at param t. If zero length, tries to move a *
  1057. * tad.                                                                       *
  1058. *                                                                            *
  1059. * PARAMETERS:                                                                *
  1060. *   Vec:     Where to place the result.                                      *
  1061. *   Crv:     Vector field to evaluate.                                       *
  1062. *   t:       Parameter value where to evaluate the vector field.             *
  1063. *   TMin:    Minimum of vector field domain.                                 *
  1064. *   TMax:    Maximum of vector field domain.                                 *
  1065. *                                                                            *
  1066. * RETURN VALUE:                                                              *
  1067. *   CagdBType:  TRUE, if successful, FALSE otherwise.                        *
  1068. *****************************************************************************/
  1069. static CagdBType EvalTangentVector(CagdVType Vec,
  1070.                    CagdCrvStruct *Crv,
  1071.                    CagdRType t,
  1072.                    CagdRType TMin,
  1073.                    CagdRType TMax)
  1074. {
  1075.     CagdRType Len,
  1076.     *R = CagdCrvEval(Crv, t);
  1077.     
  1078.     CagdCoerceToE3(Vec, &R, -1, Crv -> PType);
  1079.     Len = PT_LENGTH(Vec);
  1080.  
  1081.     if (Len < PT_NORMALIZE_ZERO && t - NORMAL_EPSILON > TMin) {
  1082.     CagdCrvEval(Crv, t - NORMAL_EPSILON);
  1083.     CagdCoerceToE3(Vec, &R, -1, Crv -> PType);
  1084.     Len = PT_LENGTH(Vec);
  1085.     }
  1086.     if (Len < PT_NORMALIZE_ZERO && t + NORMAL_EPSILON < TMax) {
  1087.     CagdCrvEval(Crv, t + NORMAL_EPSILON);
  1088.     CagdCoerceToE3(Vec, &R, -1, Crv -> PType);
  1089.     Len = PT_LENGTH(Vec);
  1090.     }
  1091.  
  1092.     if (Len > PT_NORMALIZE_ZERO) {
  1093.     Len = 1.0 / Len;
  1094.     PT_SCALE(Vec, Len);
  1095.     return TRUE;
  1096.     }
  1097.     else
  1098.         return FALSE;
  1099. }
  1100.  
  1101. /*****************************************************************************
  1102. * DESCRIPTION:                                                               M
  1103. * Evaluates the unit normals of a surface at a mesh defined by subdividing   M
  1104. * the parametric space into a grid of size UFineNess by VFineNess.         M
  1105. *   The normals are saved in a linear CagdVecStruct vector which is          M
  1106. * allocated dynamically. Data is saved u inc. first.                 M
  1107. *   This routine is much faster than evaluating normal for each point,       M
  1108. * individually.                                                              M
  1109. *                                                                            *
  1110. * PARAMETERS:                                                                M
  1111. *   Srf:          To compute normals on a grid of its parametric domain.     M
  1112. *   UFineNess:    U Fineness of imposed grid on Srf's parametric domain.     M
  1113. *   VFineNess:    V Fineness of imposed grid on Srf's parametric domain.     M
  1114. *                                                                            *
  1115. * RETURN VALUE:                                                              M
  1116. *   CagdVecStruct *:   An vector of unit normals (u increments first).       M
  1117. *                                                                            *
  1118. * KEYWORDS:                                                                  M
  1119. *   BspSrfMeshNormalsSymb, normal                                            M
  1120. *****************************************************************************/
  1121. CagdVecStruct *BspSrfMeshNormalsSymb(CagdSrfStruct *Srf,
  1122.                      int UFineNess,
  1123.                      int VFineNess)
  1124. {
  1125.     int i, j;
  1126.     CagdRType U, V, UMin, UMax, VMin, VMax, **Points;
  1127.     CagdVecStruct *Normals, *NPtr;
  1128.     CagdSrfStruct *TstNormalSrf,
  1129.     *NormalSrf = SymbSrfNormalSrf(Srf);
  1130.  
  1131.     /* Verify we have a valid regular surface. */
  1132.     TstNormalSrf = CagdCoerceSrfTo(NormalSrf, CAGD_PT_E3_TYPE);
  1133.     for (Points = TstNormalSrf -> Points, i = 0;
  1134.      i < TstNormalSrf -> ULength * TstNormalSrf -> VLength;
  1135.      i++)
  1136.     if (!APX_EQ(Points[1][i], 0.0 ) ||
  1137.         !APX_EQ(Points[2][i], 0.0 ) ||
  1138.         !APX_EQ(Points[3][i], 0.0 ))
  1139.         break;
  1140.     CagdSrfFree(TstNormalSrf);
  1141.     if (i >= TstNormalSrf -> ULength * TstNormalSrf -> VLength) {
  1142.     /* Not a regular surface - ignore it, */
  1143.     return NULL;
  1144.     }
  1145.  
  1146.     UFineNess = MAX(2, UFineNess);
  1147.     VFineNess = MAX(2, VFineNess);
  1148.  
  1149.     Normals = CagdVecArrayNew(UFineNess * VFineNess);
  1150.  
  1151.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  1152.  
  1153.     NPtr = Normals;
  1154.     for (i = 0; i < UFineNess; i++) {
  1155.     for (j = 0; j < VFineNess; j++) {
  1156.         U = UMin + (UMax - UMin) * i / (UFineNess - 1);
  1157.         V = VMin + (VMax - VMin) * j / (VFineNess - 1);
  1158.         CagdEvaluateSurfaceVecField(NPtr -> Vec, NormalSrf, U, V);
  1159.  
  1160.         NPtr++;
  1161.     }
  1162.     }
  1163.  
  1164.     CagdSrfFree(NormalSrf);
  1165.  
  1166.     return Normals;
  1167. }
  1168.  
  1169. /*****************************************************************************
  1170. * DESCRIPTION:                                                               M
  1171. * Converts a Bspline surface into a Bspline surface with floating end        M
  1172. * conditions.                                                                M
  1173. *                                                                            *
  1174. * PARAMETERS:                                                                M
  1175. *   Srf:       Bspline surface to convert to floating end conditions. Assume M
  1176. *              Crv is either periodic or has floating end condition.         M
  1177. *                                                                            *
  1178. * RETURN VALUE:                                                              M
  1179. *   CagdSrfStruct *:  A Bspline surface with floating end conditions,        M
  1180. *                     representing the same geometry as Srf.                 M
  1181. *                                                                            *
  1182. * KEYWORDS:                                                                  M
  1183. *   CnvrtPeriodic2FloatSrf, conversion                                       M
  1184. *****************************************************************************/
  1185. CagdSrfStruct *CnvrtPeriodic2FloatSrf(CagdSrfStruct *Srf)
  1186. {
  1187.     int i, j,
  1188.     UOrder = Srf -> UOrder,
  1189.     VOrder = Srf -> VOrder,
  1190.     ULength = Srf -> ULength,
  1191.     VLength = Srf -> VLength,
  1192.     MaxAxis = CAGD_NUM_OF_PT_COORD(Srf -> PType);
  1193.     CagdSrfStruct *NewSrf;
  1194.  
  1195.     if (!CAGD_IS_UPERIODIC_SRF(Srf) && !CAGD_IS_VPERIODIC_SRF(Srf)) {
  1196.     CAGD_FATAL_ERROR(CAGD_ERR_PERIODIC_EXPECTED);
  1197.     return NULL;
  1198.     }
  1199.  
  1200.     NewSrf = BspSrfNew(CAGD_SRF_UPT_LST_LEN(Srf), CAGD_SRF_VPT_LST_LEN(Srf),
  1201.                UOrder, VOrder, Srf -> PType);
  1202.  
  1203.     NewSrf -> UKnotVector = BspKnotCopy(Srf -> UKnotVector,
  1204.                     CAGD_SRF_UPT_LST_LEN(Srf) + UOrder);
  1205.     NewSrf -> VKnotVector = BspKnotCopy(Srf -> VKnotVector,
  1206.                     CAGD_SRF_VPT_LST_LEN(Srf) + UOrder);
  1207.  
  1208.     for (i = !CAGD_IS_RATIONAL_PT(Srf -> PType); i <= MaxAxis; i++) {
  1209.     CagdRType *NewPts,
  1210.         *Pts = Srf -> Points[i];
  1211.  
  1212.     NewSrf -> Points[i] = (CagdRType *) IritMalloc(sizeof(CagdRType) *
  1213.                            CAGD_SRF_UPT_LST_LEN(Srf) *
  1214.                            CAGD_SRF_VPT_LST_LEN(Srf));
  1215.     NewPts = NewSrf -> Points[i];
  1216.  
  1217.     for (j = 0; j < VLength; j++, Pts += ULength) {
  1218.         CAGD_GEN_COPY(NewPts, Pts, sizeof(CagdRType) * ULength);
  1219.         NewPts += ULength;
  1220.         if (Srf -> UPeriodic) {
  1221.         CAGD_GEN_COPY(NewPts, Pts,
  1222.                   sizeof(CagdRType) * (UOrder - 1));
  1223.         NewPts += UOrder - 1;
  1224.         }
  1225.     }
  1226.     if (Srf -> VPeriodic) {
  1227.         CAGD_GEN_COPY(NewPts, NewSrf -> Points[i],
  1228.               sizeof(CagdRType) * (VOrder - 1) *
  1229.               CAGD_SRF_UPT_LST_LEN(Srf));
  1230.     }
  1231.     }
  1232.  
  1233.     for (i = MaxAxis + 1; i <= CAGD_MAX_PT_COORD; i++)
  1234.     NewSrf -> Points[i] = NULL;
  1235.     
  1236.     return NewSrf;
  1237. }
  1238.  
  1239. /*****************************************************************************
  1240. * DESCRIPTION:                                                               M
  1241. * Converts a Bspline surface to a Bspline surface with open end conditions.  M
  1242. *                                                                            *
  1243. * PARAMETERS:                                                                M
  1244. *   Srf:       Bspline surface to convert to open end conditions.            M
  1245. *                                                                            *
  1246. * RETURN VALUE:                                                              M
  1247. *   CagdSrfStruct *:  A Bspline surface with open end conditions,         M
  1248. *                     representing the same geometry as Srf.                 M
  1249. *                                                                            *
  1250. * KEYWORDS:                                                                  M
  1251. *   CnvrtFloat2OpenSrf, conversion                                           M
  1252. *****************************************************************************/
  1253. CagdSrfStruct *CnvrtFloat2OpenSrf(CagdSrfStruct *Srf)
  1254. {
  1255.     CagdRType UMin, UMax, VMin, VMax;
  1256.     CagdSrfStruct *TSrf1, *TSrf2;
  1257.  
  1258.     if (!CAGD_IS_BSPLINE_SRF(Srf)) {
  1259.     CAGD_FATAL_ERROR(CAGD_ERR_BSP_SRF_EXPECT);
  1260.     return NULL;
  1261.     }
  1262.  
  1263.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  1264.  
  1265.     TSrf1 = CagdSrfRegionFromSrf(Srf, UMin, UMax, CAGD_CONST_U_DIR);
  1266.     TSrf2 = CagdSrfRegionFromSrf(TSrf1, VMin, VMax, CAGD_CONST_V_DIR);
  1267.     CagdSrfFree(TSrf1);
  1268.     return TSrf2;
  1269. }
  1270.  
  1271.